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Abstract 

In this paper, we present a Bayesian approach for spectral unmixing of multispectral Lidar 
(MSL) data associated with surface reflection from targeted surfaces composed of several known 
materials. The problem addressed is the estimation of the positions and area distribution of each 
material. In the Bayesian framework, appropriate prior distributions are assigned to the unknown 
model parameters and a Markov chain Monte Carlo method is used to sample the resulting posterior 
distribution. The performance of the proposed algorithm is evaluated using synthetic MSL signals, 
for which single and multi-layered models are derived. To evaluate the expected estimation per¬ 
formance associated with MSL signal analysis, a Cramer-Rao lower bound associated with model 
considered is also derived, and compared with the experimental data. Both the theoretical lower 
bound and the experimental analysis will be of primary assistance in future instrument design. 

Index Terms 

Remote sensing, Multispectral Lidar, Spectral unmixing. Estimation performance, Bayesian 
estimation, Markov Chain Monte Carlo. 

I. Introduction 

Laser altimetry (or Lidar) is an acknowledged tool for extracting spatial structures from 
three-dimensional (3D) scenes, including forest canopies 0. 0- Using time-of-flight to 
create a distance profile, signal analysis can recover tree and canopy heights, leaf area indices 
(LAIs) and ground slope by analyzing the reflected photons from a target. Conversely, passive 
multispectral (MSI) (dozen of wavelengths) and hyperspectral images (HSI) (hundreds of 
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wavelengths) are widely used to extract spectral information about the scene which, for 
forest monitoring, can also provide useful parameters about the canopy composition and 
health (tree species, leaf chlorophyll content, water content, stress, among others) 0, 0. 
The most natural evolution to extract spatial and spectral information from sensed scenes is 
to couple Lidar data and multi/hyperspectral images |5], 0. Although the fusion of Lidar 
data and HSIs can improve scene characterization, there are problems of data synchronization 
in space (alignment, resolution) and time (dynamic scene, change of observation conditions, 
etc). For these reasons, multispectral Lidar (MSL) has recently received attention from the 
remote sensing community for its ability to extract both structural and spectral information 
from 3D scenes Q,@. 

The key advantage of MSL is the ability to provide information on the vertical distribution 
of spectra, used to infer physiological processes directly linked to actual carbon sequestration 
as well as carbon stocks [9J. For example, a key capability is to directly measure and classify 
ground-based shrub infestation, which is difficult with conventional HSI as the lower spectral 
response is obscured by the ’first’ signals returned from the top of the canopy. Preliminary 
trials with existing commercial LiDAR systems with three operational wavelengths have 
already taken place [10], but as we shall demonstrate here, such systems cannot yet provide 
all the necessary spectral and structural discrimination to directly measure these processes. 

Another motivation for MSL is that HSI, even when fully synchronised, can only integrate 
the spectral response along the path of each optical ray, not measure the spectral response as 
a function of distance, e.g. depth into a forest canopy. Multiple scattering effects cannot be 
neglected in some scenes with relief or containing multi-layered objects (such as trees), 
which complicates surface detection and quantification. Most existing spectral unmixing 


(SU) techniques rely on a linear mixture assumption [11 ]—[ 151 to identify the components 
(so-called endmembers) of an image and their proportions (abundances). However, it has 
been shown that the classical linear mixing model (LMM) can be inaccurate when multiple 
scattering effects occur [ J_6), 0. Although polynomial models [18], [19) (including bilinear 
models[!20[-[23]) can be adapted for “long range” multiple scattering (in contrast to intimate 


mixtures models [24]), motivating the additional parameter constraints and designing accurate 
nonlinear unmixing methods are still challenging problems [25]. Since MSL data present an 
additional dimension when compared to HSIs, we can expect a reduction of the SU problem 
complexity and thus an improvement in the target characterization performance. 


In [26], the authors proposed a Bayesian algorithm to estimate the positions and amplitudes 
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in Lidar signals associated with a multi-layer target. This method has been extended in [;9|, 


127] to MSL by first estimating the positions of the peaks ( i.e ., the layers), which were 
assumed to be the same in all spectral bands, then estimating the amplitudes of the peaks 
for each wavelength, and finally relating the estimated amplitude peaks to areas associated 
with each material that make up the target based on a linear mixing model. The method has 
been applied to real MSL signals (four wavelengths) to analyze a conifer structure and has 
been shown to aid the recovery of bark and needle areas of the tree assuming it is modeled 
as a set of irregularly spaced layers. 

In this paper, we propose to investigate a SU problem applied to single- and multi-layered 
targets to estimate the positions and proportions of each (known) material. This problem 
extends the supervised SU problem of HSIs by also estimating the target position. Estimating 
the spectral variation of material signatures (unsupervised SU problem) is also of interest. 
For example, the estimation of physiological content, such as Chlorophyll concentration, in 
forest canopies is of key concern, and this results in variable spectra for leaves and needles. 
Extracting these spectra from MSL signals is a more challenging problem (probably more 
difficult than in HSIs due to the statistical properties of the noise and observation model) 
which is out of scope of this paper, but will require particular attention in future work. 

In contrast with the classical additive Gaussian noise assumption used for HSIs, a Poisson 
noise model is more appropriate for MSL signals. Indeed, Lidar and thus MSL systems usually 
record, for each pixel/region of the scene, a histogram of time delays between emitted laser 
pulses and the detected photon arrivals. Within each histogram bin, the number of detected 
photons follows a discrete distribution which can be approximated by a Poisson distribution 
due to the particle nature of light. Using a Bayesian approach, appropriate prior distributions 
are chosen for the unknown parameters of the model considered here, i.e., the material areas, 
the surface positions and the background parameters. The joint posterior distribution of these 
parameters is then derived. Since the classical Bayesian estimators cannot be easily computed 
from this joint posterior (mainly due to the model complexity), a Markov chain Monte Carlo 
(MCMC) method is used to generate samples according to the posterior of interest. 

More precisely, following the principles of the Gibbs sampler, samples are generated 
according to the conditional distributions of the posterior. Due to the possibly high correlations 
between the material proportions/areas, we propose to use a Hamiltonian Monte Carlo (HMC) 


128] method to sample according to some of the conditional distributions. HMCs are powerful 
simulation strategies based on Hamiltonian dynamics which can improve the convergence 
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and mixing properties of classical MCMC methods (such as the Gibbs sampler and the 


Metropolis-Hastings algorithm) [29], [301. These methods have received growing interest in 


many applications, especially when the number of parameters to be estimated is large [3TJ, 
p2) . Classical HMC can only be used for unconstrained variables. However, new HMC 
methods have been recently proposed to handle constrained variables 


Chap. 5], [33), [34) 

which allow HMCs to sample according to the posterior of the Bayesian model proposed for 
SU. Finally, as in any MCMC method, the generated samples are used to compute minimum 
mean square error (MMSE) estimators as well as measures of uncertainties such as confidence 
intervals. 

Predicting the parameter estimation performance is of prime interest for designing an 
estimation procedure but also assists with the instrument design. Indeed, since MSL is a recent 
modality, it is important to identify the key parameters that have an influence on the material 
estimation performance (such as the necessary signal to background levels, the number of 
bands to be considered and the associated wavelengths). To guide future instrument design, 
we consider a Cramer-Rao lower bound (CRLB) associated with the observation model and 
show that it can be used to approximate the estimation errors of the proposed Bayesian 
algorithm. 

The remainder of the paper is organized as follows. Section [IT] introduces the observation 


model associated with MSL returns for a single-layered object to be analyzed. Section III 


presents the hierarchical Bayesian model associated with the proposed model and its posterior 
distribution. The hybrid MCMC method used to sample from the posterior of interest is 
detailed in Section [IVj Section [V] investigates the expected SU performance using Cramer- 


Rao lower bounds. Experimental results are shown and discussed in Section VI Conclusions 
are reported in Section |VII 


II. Problem formulation 

This section introduces the observation statistical model associated with MSL returns 


for a single-layered object which will be used in Sections III and IV to solve the SU 
problem of MSL data. Precisely, we consider a set of L observed Lidar waveforms = 
[ye t i,..., ye,T) T , £ £ {1,..., L} where T is the number of temporal (corresponding to range) 
bins. To be precise, y^ t is the photon count within the /1h bin of the (th spectral band 
considered. Let t 0 be the position of a complex object surface or surfaces at a given range from 
the sensor, composed of R materials whose spectral signatures (observed at L wavelengths) 
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are denoted as m r = [m l r ,..., m ^ r } T , r e {1,..., f?}. According to [26|, each photon count 
yvt is assumed to be drawn from the following Poisson distribution 


Ue,t ~ 'P 


R 

^ ^ ^r‘^T'£,r90,c(j' S ) 

r =1 



(i) 


where go A') is the photon impulse response whose shape can differ between wavelength 
channels, w r denotes the area of the rth material composing the object (relative to a reference 
area) and be is a background and dark photon level, which is constant in all bins at a given 
wavelength. Without loss of generality the photon impulse responses are assumed to be the 
same and modeled by the following piece-wise exponential 


' 

go At) = /3 < 




_ T 1 t + Tl-tg 

exp 2^ exp T i if t — t 0 < — Ti 

tt-tp > 2 

exp 2^2 if T\ < t — t o < Ti 

_ t 2 t-T-2-tQ 

exp 2^ exp T 2 if T 2 < t — t 0 < T 3 

_ t 2 Tz-t* 

exp 2?^ exp T 2 exp T 3 else 


where <fi = [Xj, T 2 , T 3 , Ti, r 2 , r 3 , a 2 , (3] T is a positive hyperparameter vector. In this study, 
we assume that the shape parameters of the model (which is appropriate for our own lidar 
sensors) are fixed and known from the instrumental response. As seen later in Fig. [3] the shape 
is slightly asymmetrical when compared to the more common Gaussian model [35] which 
is used for higher power, longer pulse duration lidar systems. The consideration of a piece- 
wise exponential approximation is motivated by the actual shape of our instrumental impulse 
response but it is worth mentioning that any positive analytical approximation, even different 
approximations for the different wavelengths, could be used instead without modifying the 
proposed Bayesian estimation procedure. The consideration of a single impulse response 
model is done only for ease of reading. While a Gaussian approximation can lead to a slight 
bias in depth estimation, we shall use a Gaussian approximation in Section [V] to compute 
the Cramer-Rao bounds since the Gaussian approximation is twice differentiable (in contrast 
to the piece-wise exponential approximation), which simplifies the derivation of the bounds. 

Due to physical considerations, the relative areas are assumed to satisfy the following 
positivity constraints 


w r > 0, r = 1,..., R. 


( 2 ) 


Similarly, the average background level in each channel satisfies be > 0, i — 1,..., L. Spectral 
unmixing of the MSL data consists of estimating the position ( i.e ., t 0 ) of the target, the area 
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w r of each component m r of the target, as well as the background levels from the observed 
data in Y = [yi,..., yj. The next section studies a Bayesian model to estimate the unknown 
parameters in ([I]) while ensuring the positivity constraints mentioned above. 

III. Bayesian model 

The unknown parameter vector associated with (jT]) contains the relative areas w = [w i...., wr] t , 
the position of the target t 0 and background levels b = [ 61 ,... ,b L ] T (satisfying positivity 
constraints). This section summarizes the likelihood and the parameter priors associated with 
these parameters. 

A. Likelihood 

Assuming that the detected photon counts/noise realizations, conditioned on their mean in 
all bins and for all wavelengths, are independent, Eq. ([[]) leads to 

L T yVl,t 

P(Y|M, w,b,f 0 ) — Till ~^ 7 ex P Af,< (3) 

where M = [nix,..., m R ], \^ t = m^.yfg Q ^{t) + be and m^. denotes the fth row of M. 

B. Prior for the target position 

Since we don’t have prior information about the position of the target, the following 
uniform distribution 

to ~ W(l;T)(fo) (4) 

is assigned to t 0 . Note that the position to is a real variable that is not restricted to the integer 
values in (1; T). 

C. Prior for the relative areas 

To reflect the lack of prior knowledge about the areas, the following truncated Gaussian 
prior 

w r \a 2 ~ A/k+(uv; 0, a 2 ) (5) 

is assigned to each area w r , where A/r+(-; 0, a 2 ) denotes the Gaussian distribution restricted 
to M + , which hidden mean and variance 0 and a 2 , respectively. The hyperparameter a 2 is 
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shared by all the parameters w r and is arbitrarily fixed to a large value to ensure a weakly 
informative prior. Assuming prior independence for the unknown areas yields 

f{w\a 2 )(xf-^j exp~^ wTw l (R+) fl (w) (6) 

where oc means ’’proportional to” and (•) is the indicator function defined on (E + ) /l> . 

D. Prior for the background level 

Similarly, assigning a truncated Gaussian prior to each background level and assuming 
prior independence between these parameters leads to 

L 

/( b l7 2 ) °c exp“^ b h l {R+)L (b) (7) 

where the hyperparameter 7 2 is shared by all the parameters hi and is arbitrarily fixed to a 
large value to ensure a weakly informative prior. 

In this paper, we assume that prior knowledge about the relative areas and background 
levels is limited. We use weakly informative Gaussian priors restricted to the positive orthant 
to reflect this lack of knowledge and reduce the estimation bias while ensuring the positivity 
constraints are satisfied. However any proper weakly informative prior distributions could 
have been used, leading to similar estimation results. Similarly, if useful information about 
the unknown parameters is available, more informative (possibly physics-based) priors can 
be used instead of Q and ((6)) (target position and composition) and ([7j) (background levels). 

E. Joint posterior distribution 



Fig. 1. DAG for the parameter and hyperparameter priors (the fixed parameters appear in boxes). 


The resulting directed acyclic graph (DAG) associated with the proposed Bayesian model 
is depicted in Fig. [Tj The joint posterior distribution of the unknown parameter vector 6 = 
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{w,b,t 0 } can be computed using 


f(0 1Y, M, a, P) oc P(Y|M, w, b,t o )f(0\a, (3) (8) 

where P(Y|M, w,b,£ 0 ) has been defined in ([3]) and f(0\a,(3) = f(w\a 2 )f(b\/3 2 )f(t 0 ). 

Unfortunately, it is difficult to obtain closed form expressions for standard Bayesian esti¬ 
mators associated with ([8]). In this paper, we propose to use efficient Markov Chain Monte 
Carlo (MCMC) methods to generate samples asymptotically distributed according to ([8]), 
employing a Gibbs sampler. The principle of the Gibbs sampler is to sample according to 
the conditional distributions of the posterior of interest J30[ Chap. 10]. Due to the high 
correlation between the elements of w (which will be discussed later in the paper), we use 
an Hamiltonian Monte Carlo (HMC) method which improves the mixing properties of the 
sampling procedure (compared to a classical Metropolis-within-Gibbs sampler using Gaussian 
random walks). This method is detailed in the next section. 

IV. Bayesian inference using a Constrained Hamiltonian Monte Carlo 

method 

In this Section, we develop an HMC-within Gibbs sampler to generate samples according 
to ([8]) and estimate the unknown parameters involved in ([[} in order to solve the SU problem 
of MSL data for a single-layered object. The proposed sampler consists of three steps to 
update sequentially w, t 0 and b and is summarised in Algo. [T} 

A. Sampling the areas 

The full conditional distribution of w is given by 

/(w|Y, M, t 0 , b, a 2 ) oc P(Y|M, w, b, £ 0 )/(w|a 2 ) (9) 

and is not a standard distribution which is easy to sample. Consequently, it is the norm to 
use an accept/reject procedure to update b. A classical and simple approach would use a 
multivariate Gaussian random walk. However, in practice the elements of w can be highly 
correlated, especially when the materials present similar spectral signatures (which will be 
the case for vegetation targets). Consequently, a Hamiltonian Monte Carlo method [ |29| Chap. 
5] is preferred to improve the mixing properties of the sampler. The principle of HMCs is 
to introduce auxiliary, or momentum, variables, and perform a Metropolis-Hastings move in 
a higher dimensional parameter space. The proposal distribution is then built to take into 
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account the shape of the target distribution ([9]). Due to the area constraints Q, a constrained 
HMC must be used. In this paper, we used a constrained HMC (CHMC) scheme similar to 


that described in [19] and thus consider the following potential energy function 

U{ w) = - y e , t log(A^) + A e , t + 


T 

w w 


Lt 


2a 2 


( 10 ) 


to simulate Hamiltonian dynamics and compute the appropriate acceptance ratio (see [19|, 
[29] for technical details). This choice of potential energy function ensures that U{ w) = 
— log (/(w|Y, M, t 0 , b, a 2 )) + c where c is a positive constant. Note that the proposed CHMC 
step can be applied for any differentiable relative area prior. Metropolis-Hastings or more 
complex HMC moves should be used instead when considering non-differentiable priors 
instead of ([6]). 


B. Sampling the target position 
It can be shown from ([8]) that 

f(to\Y, M. w, b, a 2 ) oc P(Y|M. w, b, t Q )f(t 0 ) 


( 11 ) 


is not a standard distribution and an accept/reject procedure must be used to update the 
target position t 0 . We use a Gaussian random walk to update this parameter. More precisely, 
a doubly truncated Gaussian proposal is considered to ensure each candidate belongs to the 
admissible set (1; T ) and the variance of the proposal is adjusted during the burn-in period 


of the sampler to obtain an acceptance rate close to 0.45, as recommended in [36, p. 8]. 


C. Sampling the background levels 
It can be shown from ([8]) that 

L 

/(b|Y, M, w, t 0 , « 2 ) = f{be\ye, w, t 0 , a 2 ), (12) 

£= 1 

where 

f(be\yei m^ >: , w, t 0 , a 2 ) 

T b 2 

°C Aexp _A£ ’* exp - ^ 1 M+ (b e ) (13) 

t= l 

Sampling from ( [13] ) is again not straightforward and Gaussian random walks are used to 
update the background levels, similar to the position t 0 . However the background levels are 
a posteriori independent and can be updated in parallel. Similar to the target position update, 
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the variances of the L parallel Gaussian random walk procedures are set during the burn-in 
period of the sampler to obtain an acceptance rate close to 0.45. 

Algorithm 1 

Gibbs Sampling Algorithm (single layer) 

1: Fixed input parameters: M, </>,a 2 , 7 2 , number of burn-in iterations N hl , total number of iterations)V M c 
2: Initialization (i = 0) 

• Set w(°), t£j°\ b® 

3: Iterations (1 < i < Nmc) 

4: Sample from the pdf in (|9]i and CHMC 
5: Sample from the pdf in CD and a Gaussian random walk 
6: Sample bW from the pdfs in ( f]~3j ) and Gaussian random walks 
7: Set i = i + 1. 

After generating A r MC samples using the procedures detailed above and removing N bi 
iterations associated with the burn-in period of the sampler (N b \ has been set from preliminary 
runs), the MMSE estimator of the unknown parameters can be approximated by computing 
the empirical averages of the remaining samples. The minimal length of the burn-in period can 


be determined using several convergence diagnostics [36 ] but the number of initial samples to 
be discarded generally varies between data sets and can be thus difficult to adjust in advance 
without overestimating it. In previous work on parallel acceleration of RJMCMC algorithms 
we have investigated the application of such diagnostics to our monochromatic LiDAR signals 
[371, but as absolute speed of processing is not a major concern here, the length of the burn-in 
period is assessed visually from the preliminary runs and fixed to ensure that the sampler has 
converged. The total number of iterations has then been set to obtain accurate approximations 


(computed over 4000 samples in Section VI) of the MMSE estimators. 


V. Predicting unmixing performance 

This section studies a Cramer-Rao lower bound associated with the observation model 
(|Tj) which can be used to assess the performance of methods that aim to solve the SU 
problem of MSL data considered in this paper as well as to assist with the future instrument 


design. Precisely, Section V-A recalls the definition of the CRLB for the problem of interest 


and Section V-B discusses the impact of several key parameters on the expected parameter 


estimation performance of SU methods. 
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A. Cramer-Rao lower bound 

Prediction of the unmixing performance is necessary to assist in the design of a lidar 
system for a specific application, identifying those parameters that have the most significant 
impact. In deriving a CRLB, we propose firstly to relax the impulse response model in ([2]) 
by considering the following Gaussian approximation 

(t-t 0 ) 2 

goA t )=P ex P 2ct2 Vt G ( 1 ; T), W. (14) 


This simplifies the CRLB derivation and does not significantly bias the prediction as the piece- 
wise exponential impulse response can be accurately approximated by a Gaussian function. 
The CRLB associated with any unbiased estimator 6 of the parameter vector 6 involved in 
the mixing model (|T|) (having replaced by 9o/(t)) and constructed from Y is given by 

CRLB(0) = Jp 1 

where J F is the Fisher information matrix whose elements areQ 

'd 2 log /(Y|0)" 


(15) 


“ _E Y|0 


d9.;d6i 


ihJ — 1, • • •,R + L + 1. 


The ith diagonal element of the CRLB matrix in ( fl5l ) provides a lower bound for the 
variance of 9 t , given that 0 is an unbiased estimator of 6. Of course, the Bayesian estimation 


procedure proposed in Sections III and IV does not provide a strictly unbiased estimator 
of 0 and a Bayesian CRLB should have been considered instead of (p3|). However, due to 


the weakly informative priors chosen in Section III, when the actual vector 6 is far enough 
from the boundaries of the admissible set (defined by the positivity constraints), the bias 
of the proposed estimation procedure can be neglected and the resulting estimator achieves 
the CRLB (as will be shown in Section |VT|). Moreover, the CRLB in ( [15] ) can also provide 
information about the estimation performance of a potential optimization method that could 
be proposed to estimate 0 based on maximum likelihood estimation (MLE). 


B. Performance analysis 

We first investigate the performance of spectral unmixing of MSL signals for a single 
layered, artificial target composed of R — 3 materials, specifically needles, bark and soil. 
These materials have been chosen because of our interests in forest canopy monitoring using 
MSL signals [9). The reflectance spectra of these materials, observed at equally spaced 

'The Fisher information matrix Jf is derived in the Appendix 
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Fig. 2. Spectral signatures of needle (red), bark (blue) and soil (green) considered in the synthetic data. 


spectral bands ranging from 400nm to 2500nm are depicted in Fig. [2j It is interesting to note 
the strong similarity between the bark and needles spectral signatures, which can make their 
discrimination difficult and highlights the need of multiple wavelengths to quantify these 
materials. The maximum number of spectral bands considered in this paper has been set to 
L max = 32, which is a realistic value for a short-term real measurement campaign (the current 
instrument uses only 4 wavelengths). The relative area of needles (resp. bark and soil) has 
been arbitrarily set to w n = 0.2 (resp. w b = 0.3 and w s = 0.4). These relative areas do not 
necessarily sum to one as we consider that part of the laser light can penetrate through the 
target (semi-transparent target or a significant gap fraction) and may not be further reflected 
(single hit assumption). The fixed model parameters have been fixed from experiments to 


( 16 ) 


T = 2500 
P = 3000 
to = 1000 
b t = b = 10, W 

These parameters will be fixed in the remainder of the paper unless otherwise specified. The 
impulse response parameters have been set to T\ = 402, T 2 = 12.5, T 3 = 239, T\ = 395, r 2 = 
7.9, r 3 = 1595 and cr 2 = 105.82 for the piece-wise exponential approximation in ([2]) and 
cr 2 = 105.68 for the Gaussian approximation in ( [T4| ). These parameters have been obtained 
by fitting the experimental impulse response measured in [261. Fig. [3] shows that the Gaussian 
approximation provides a good estimate of the experimental impulse response, although the 
piece-wise exponential better fits the experimental curve. Fig. [4] shows an example of MSL 
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Fig. 3. Real normalized photon impulse response (red dots) and approximations using piece-wise exponential (blue line) 
and Gaussian (black line) functions. 


data generated using the parameters in ( [T6| ) and the impulse response approximation in Q. 
Table |T] shows the predicted estimation performance for the unknown parameters of interest 
(; i.e ., w = [w n ,w b ,w s ] T and t 0 ) assuming the Gaussian approximation of the photon impulse 
response ©. The relative errors provided in the bottom row of Table |T] are computed by 
dividing the square root of the CRLBs by the actual values of the parameters. The relative 
estimation errors of the background levels (not presented here) are lower than 1%. 


TABLE i 

Estimation performance (L — 32). 



Wn 

w b 

W s 

to 

Actual value 

0.2 

0.3 

0.4 

1500 

CRLB 

2.6 x 10" 4 

0.001 

3.6 x 10" 5 

1.8 x 10“ 4 

Rel. err. (%) 

8.0 

10.7 

1.5 

< 10" 5 


These results show that the estimation errors associated with the target position are usually 
much lower than those associated with the material areas. Moreover, although the amplitude 
of the peak for each wavelength is much larger than the background level (see Fig. [4]), the 
CRLB predicts possibly large estimation errors, especially for the needle and bark areas. 
This can be partially explained by the fact that the soil reflectance spectrum presents an 
average energy which is higher than those of needles and bark (and thus w s > Wb > w n ) 
but also and mainly because of the high correlation between the bark and needle spectra, 
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Fig. 4. Example of MSL signals for L = 32 spectral bands and a target located at to = 1000 (the bottom subplot is a 
zoom around the actual target location). The different colors correspond to different spectral bands. 


which complicates their quantification. Figs. [5] to [7] show the predicted estimation errors 
for different values of the key parameters /3 (which is related to the number of photons 
emitted by the laser sources), the number of spectral bands (equally spaced between 400 nm 
and 2500 nm) and the average background level b (assumed to be the same in all spectral 
bands and bins). Fig. [5] shows that the estimation performance generally increases with the 
number of spectral bands. This result is well known for SU of multispectral and hyperspectral 
images and is here demonstrated for MSL signals. However, the performance improvement 
can vary, depending on the additional spectral bands considered. For instance, additional 
bands containing similar spectral information may not significantly improve the unmixing 
results ( e.g ., several wavelengths between 2000 nm and 2500 nm for the three materials in 
Fig. |2j. Fig. [6] shows that the parameter (3 has a significant impact on the SU performance. 
This parameter increases with the amplitude and the number of laser source pulses used 


October 6, 2015 


DRAFT 



















15 


to acquire the data and decreases with the average distance between the source and the 
target (as the probability of recording a reflected photon decreases). Increasing f3 leads to 
better area estimation but can require a longer target exposure (which can be problematic 
for airborne sensors for instance). Finally, Fig. [7] shows that for relatively small values of 
background levels (compared to (3), the SU performance is not significantly degraded when 
the background increases. This background level (which depends primarily on the background 
(e.g. solar) radiation as well as the instrument design) is expected to be quite small in practice. 



Fig. 5. Evolution of the CRLB for w and to as a function of the number of bands L (/3 = 3000, bi = 10, W). 



Fig. 6. Evolution of the CRLB for w and to as a function of /3 (L = 32, bt = 10). 
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Fig. 7. Evolution of the CRLB for w and to as a function of b = be,\/l (L = 32, /? = 3000). 
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VI. Experiments 

In this section, we first apply (Section |VI-A[ ) the proposed SU method to synthetic MSL 
data associated with a single-layered target and compare its estimation performance to that 
predicted by the CLRB presented in Section [V] and to that of a state-of-the art method [9]. 


In Section VI-B we investigate the case of a multi-layered target whose layer positions are 
assumed to be known. Although the current algorithm does not estimate the positions of 


multiple surfaces, Section VI-B aims to illustrate how the proposed method can be extended 
to analyze more complex objects. The main steps of the methodology (proposed to analyze 
single and multi-layer target) are summarised in Fig. [8j 



Fig. 8. Main steps of the method proposed for spectral unmixing of MSL data associated with single and multi-layer 
targets. 


A. Single layer target 

In this section, we evaluate the performance of the proposed spectral unmixing algorithm 
on synthetic MSL signals generated using the parameters used in Section |V-B| (Eq. < [T6| )) and 
impulse response approximation Q. 

The number of spectral bands has been chosen as L = 4,8,16, 32, equally spaced between 
400 nm and 2500 nm. For each scenario, the number of iterations has been set to N MC = 8000 


October 6, 2015 


DRAFT 


















18 


(including N bi = 4000 bum-in samples). The hyperparameters have been fixed to (« 2 ,7 2 ) = 
(10 6 ,10 6 ). The estimation performance of the proposed algorithm is evaluated by comparing 

the CRLB defined in Section [V] to the mean square error (MSE) defined as 

' / A \ 2" 


MSE, 


0i — Qi 


i — 1,..., K T f f 1 


(17) 


where 6^ and 0, are the actual and estimated /1h element of the unknown parameter vector 
6 and the expectation is approximated using 900 Monte Carlo runs. The performance of 
the proposed algorithm is compared to the CRLB (including the Gaussian approximation 
( fT4| )) and to the performance of the algorithm developed in [9j, denoted as “State-of-the-art” 
and which consists of estimating sequentially the position, amplitudes of the peaks for each 
wavelength and finally the relative areas. 

Table [11] shows that the MSEs obtained for L 6 {4,8,16,32} by the two algorithms are 
close to the CLRBs, although the proposed algorithm generally outperforms the method 
proposed in [0 due to the joint estimation of the target position and the material areas. Note 
that for L — 4, the variance of the proposed estimator is slightly lower than the CRLB. This 
can be mainly explained by the fact that the estimator is no longer unbiased in that case. 


Table III compares the average relative estimation errors obtained by the two algorithms and 
computed by dividing the square root of the MSEs by the actual values of the parameters. 
These results show that increasing the number of spectral bands from L = 4 to L = 32 
almost divides the area estimation errors by three (e.g., from « 30% to « 10% for the bark 
area), which highlights the benefits of increasing the number of wavelengths. 


TABLE II 

Estimation performance: MSE . 



w n 

w b 

w s 

Actual parameter value 

0.2 

0.3 

0.4 

L = 32 

CRLB 

2.6 X 10“ 4 

1.0 X 10“ 3 

3.6 X 10“ 5 

Proposed Algo. 

2.7 X 10“ 4 

1.1 X 10“ 3 

3.8 X 10“ 5 

State-of-the-art 

3.0 X 10“ 4 

1.2 x 10 -3 

4.1 x 10“ 5 

L = 16 

CRLB 

4.6 X 10“ 4 

1.9 X 10 -3 

6.8 X 10 -5 

Proposed Algo. 

4.8 x 10“ 4 

2.0 X 10“ 3 

7.0 X 10 -6 

State-of-the-art 

5.3 X 10“ 4 

2.2 x 10 -3 

7.8 X 10 -6 

II 

00 

CRLB 

5.8 X 10“ 4 

2.6 X 10 -3 

1.0 X 10“ 4 

Proposed Algo. 

5.7 X 10“ 4 

2.5 X 10“ 3 

1.0 X 10“ 4 

State-of-the-art 

6.1 X 10“ 4 

2.7 X 10“ 3 

1.1 X 10“ 4 

II 

CRLB 

1.8 X 10“ 3 

8.0 X 10“ 3 

3.0 X 10“ 4 

Proposed Algo. 

1.6 X 10 -3 

7.4 X 10“ 3 

2.9 x 10“ 4 

State-of-the-art 

1.8 X 10 -3 

8.6 X 10“ 3 

3.3 X 10“ 4 
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TABLE III 

Estimation performance: Average relative errors (in %). 



w n 

w b 

w s 

Actual parameter value 

0.2 

0.3 

0.4 

L = 32 

Theo. (%) 

8.0 

10.7 

1.5 

Proposed Algo. 

8.2 

11.0 

1.5 

State-of-the-art 

8.6 

11.5 

1.6 

L = 16 

Theo. (%) 

10.8 

14.6 

2.0 

Proposed Algo. 

11.0 

14.8 

2.1 

State-of-the-art 

11.6 

15.7 

2.2 

L = 8 

Theo. (%) 

12.1 

16.8 

2.5 

Proposed Algo. 

12.0 

16.8 

2.5 

State-of-the-art 

12.4 

17.4 

2.6 

L = 4 

Theo. (%) 

21.0 

29.9 

4.3 

Proposed Algo. 

20.0 

28.8 

4.2 

State-of-the-art 

21.4 

30.8 

4.5 


B. Extension to a multi-layer target 


Algorithm 2 

Gibbs Sampling Algorithm (multi-layer) 

1: Fixed input parameters: M, </>,a 2 , 7 2 , number of layers D, layer positions {td} d=1 D , number of 

burn-in iterations TVy, total number of iterationsA\ic 

2: Initialization ( i - 0) 

. Set {) , t/°) 

1 J d=l,...,D 

3: Iterations (1 < i < Nuc) 

4: for d=l:D do 

5: Sample w rf from its conditional distribution and CHMC 

6: end for 

7: Sample bW from its conditional distribution and Gaussian random walks 
8: Set i = i + 1. 

In this section we extend the model O by considering a target composed of D layers, 
leading to 

Vex ~ V m es w d9dj{t) + be j (18) 

where gd.e(r) is the photon impulse response of the c/th layer located at t d and w () = 
[wi t d, ■ ■ ■, u'R.d] 1 h 0 is the relative area vector associated with the c/th layer. In this scenario, 
we assume that the number and positions of the layers are known. Consequently, the unmixing 
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problem reduces to estimating the relative areas of the R known components for the D 
layers (and the background levels). The joint spectral and spatial unmixing problem is more 


challenging and will be discussed in Section VI-C The Bayesian model presented in Section 
In] has been extended by considering the same prior (|6]) for the area vectors of each layer. The 

has been modified in order to update sequentially 


sampling procedure studied in Section | 
the D area vectors and the target position update step has been removed since the surface 
positions are assumed to be fixed in this paragraph (see Algo. [2]). 

1) Target description: We evaluate the unmixing performance of the extended Bayesian 
algorithm using an artificial target composed of D = 3 layers. More precisely, the multi-layer 
target is assumed to be far enough from the source to ensure that the laser rays hitting the 
target have the same incident direction. Thus the photons reflected onto the dth layers lead 
to the same impulse response ( i.e ., same tf). The first two layers are located at t\ = 1000 
and t -2 = 1500 and are composed of needles, bar and soil. It is important to note that the 3 
materials composing the first two layers are assumed to randomly distributed, without overlap. 
The third layer modeling the reference spectralon is located at f 3 = 2000. The associated 
relative areas are presented in Table [IV] and the location of the D — 3 layers is depicted in 
Fig. [9] Note that in this particular scenario, the relative areas correspond to areas visible by 
the source (no occlusion) and satisfy X^z-i w r,d = 1- However, this constraint is not 

ensured by the proposed estimation procedure. MSL data associated with this scene have 
been generated according to ( fT8| ) with (3 = 10 4 , b e = b = 10, VC and L = 32 and are depicted 
in Fig. [TO] 


2) Estimation performance: The proposed Bayesian algorithm extended to account for 
multiple layers has been applied to the MSL signals with N MC = 8000 (including A r bi = 4000 


burn-in samples). The layer positions have been fixed to their actual values. Table IV shows 
the estimated relative areas for each layer and L <E {4,8,16,32}. This table shows that the 
proposed algorithm provides accurate area estimates for each layer and the more wavelengths 


the better the estimation, as already observed with the single-layered target in Section VI-A 


C. Towards real data analysis 

Due to the current limitations of the instrument, the proposed algorithm has only been 
applied to synthetic MSL signals whose underlying model has been shown to be in good 
agreement with our previous real MSL measurements using fewer wavelengths (91, [ ]26| . 
Using this real data, we have managed to investigate the precision of area estimation in 
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Incident beams 




V 



1400 - 


1600- 


1800 - 



Fig. 9. Artificial 3-layered target composed of needles (red), bark (green), soil (blue) and pure reflective material (black). 



Fig. 10. MSL signals (L = 32) associated with the artificial 3-layered target (ti = 1000, £2 = 1500, fi = 2000). The 
different colors correspond to different spectral bands. 


a multi-layer model; ground truth was available for the endmembers, but unfortunately we 
did not have structural truth. As stated above, real MSL signal acquisition campaigns are 
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TABLE IV 

Estimation performance: Multi-layer target. 



Wn 

w b 

w s 

spectralon 

Actual 

Layer 1 

0.099 

0.099 

0.102 

0 

Layer 2 

0.080 

0.200 

0.120 

0 

Layer 3 

0 

0 

0 

0.30 

Estimated 

(£=32) 

Layer 1 

0.0928 

0.0780 

0.1185 

0.0004 

Layer 2 

0.0669 

0.1953 

0.1161 

0.0006 

Layer 3 

0.0019 

0.0011 

0.0011 

0.3274 

Estimated 

(£=16) 

Layer 1 

0.0932 

0.0773 

0.1181 

0.0006 

Layer 2 

0.0620 

0.1961 

0.1179 

0.0006 

Layer 3 

0.0024 

0.0018 

0.0014 

0.3270 

Estimated 

(£=8) 

Layer 1 

0.0957 

0.0764 

0.1166 

0.0008 

Layer 2 

0.0640 

0.1954 

0.1172 

0.0006 

Layer 3 

0.0024 

0.0028 

0.0016 

0.3268 

Estimated 

(£—4) 

Layer 1 

0.1286 

0.0687 

0.1013 

0.0015 

Layer 2 

0.0975 

0.1863 

0.1014 

0.0016 

Layer 3 

0.0023 

0.0019 

0.0015 

0.3271 


underway with 3 wavelengths [ 10] but we need to extend these and have simultaneous ground 
truth measurements to better evaluate the accuracy of the model and the estimation algorithm 
presented in this paper. Moreover, the number and positions of the layers which the multi¬ 
layer target is composed of were assumed to be known in Section VI-B] Estimating these 
parameters, especially the number of layers, for MSL data is a challenging problem that can 
be addressed in the Bayesian framework using reversible jump MCMC methods, in a similar 


fashion to [26]. This issue is however out of scope of this work and will also be addressed 
in a future paper. 


VII. Conclusion 

In this paper, we have proposed and developed a Bayesian model and an MCMC method 
for spectral unmixing of multispectral Lidar data. First, we evaluated our method on a single 
layer simulated target composed of known materials so that the MSL returns consisted of 
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the sum of the individual contributions of the different components. The algorithm estimated 
the target position, the relative areas of the materials and the noise statistical properties. 

The Cramer-Rao lower bound associated with the observation model was derived to identify 
the key parameters involved in the expected performance of the SU algorithm. It was shown 
that this bound can be used to help design future multispectral/hyperspectral Lidar systems 
(e.g., number of spectral bands and associated wavelengths, laser power ...) for specific 
applications. It was also shown that the performance of the proposed SU strategy is close 
to the results provided by the Cramer-Rao lower bound, although the bound considered in 
the paper only holds for deterministic parameters. The proposed Bayesian model was then 
extended to handle multi-layer targets (assuming the layer positions are known) and the 
simulation results provided interesting results for solving the more challenging joint spectral 
and spatial unmixing problem. 

The model use in this paper does not take into account possible multiple scattering effects 
(which are likely to occur when analyzing multi-layer scenes). Such effects have already 
been observed in HSIs and Lidar signals over canopies. Developing non-linear models for 
MSL signal analysis is a challenging problem to be addressed in future studies. 

Appendix: Fisher information matrix 

The likelihood of the observation matrix Y can be expressed as 

T L / \ vyz t 

/(Y|w,M,i 0 ,b) =TTTTEyi_exp- J « (19) 

vtf'- 

where \^ t = m^ : w go,e(t) + be- The corresponding log-likelihood P = log/(Y|w,M, t 0 ,b) 
is given by 

T L 

yt,t log (A itt ) - log (ye/) - A , l>t . (20) 

t =i t =i 

The partial derivatives of P with respect to (w.r.t.) the unknown model parameters are 

dP 
dw 

dP 
dbe 

dP 
d t 0 


L T 

EZ 

i=i t= i 
T 

ST' Vi,t _ i 

~i A e,t 


ye,t.go,e(t) 


A 


e,t 


- go/(t) m£ : 


t =i 
T L 


t= 1 1=1 


m U w Mt) 


t-tn 


17 z 
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Straightforward computations lead to 

'd 2 P' 


dw 2 

d 2 P 

dbe.dbt' 

' d 2 P 
dwdbe 

v d 2 P 


K 


d 2 P 

dt 0 dw 

d 2 P 

dt 0 dbi 


t L ~2 


Y^ 9o,e(t) T 


t =i e=i 


e,t 


«e = 


else 


9oAt) 

1 - - 


t =1 
T L 

EE 

t= 1 €=1 
T L 




[(t - t 0 )me-wg 0/ (t)Y 


cr 4 A. 


T,t 


(t - f 0 )m £ , ; w^ (t) 

- - PP -— 


t =1 £=1 
T 


E 

t=i 


cr 2 A^ ;t 

(f - f 0 )m £) :W^o^(f) 


<j 2 A 


c,t 


using the fact fact E [y^ jt ] = m^ : w^ 0 ^(^) + 
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